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Abstract — The emergence of synchronization in a network 
of coupled oscillators is a pervasive topic in various scientific 
disciplines ranging from biology, physics, and chemistry to 
social networks and engineering applications. A coupled oscilla- 
tor network is characterized by a population of heterogeneous 
oscillators and a graph describing the interaction among the 
oscillators. These two ingredients give rise to a rich dynamic 
behavior that keeps on fascinating the scientific community. 
In this article, we present a tutorial introduction to coupled 
oscillator networks, we review the vast literature on theory and 
applications, and we present a collection of different synchro- 
nization notions, conditions, and analysis approaches. We focus 
on the canonical phase oscillator models occurring in countless 
real-world synchronization phenomena, and present their rich 
phenomenology. We review a set of applications relevant to 
control scientists. We explore different approaches to phase 
and frequency synchronization, and we present a collection of 
synchronization conditions and performance estimates. For all 
results we present self-contained proofs that illustrate a sample 
of different analysis methods in a tutorial style. 

I. Introduction 

The scientific interest in synchronization of coupled oscil- 
lators can be traced back to the work by Christiaan Huygens 
on "an odd kind sympathy" between coupled pendulum 
clocks [1], and it still fascinates the scientific community 
nowadays [2], [3]. Within the rich modeling phenomenology 
on synchronization among coupled oscillators, we focus on 
the canonical model of a continuous-time limit-cycle oscil- 
lator network with continuous and bidirectional coupling. 

A network of coupled phase oscillators: A mechanical 
analog of a coupled oscillator network is the spring network 
shown in Figure [T] and consists of a group of kinematic 
particles constrained to rotate around a circle and assumed 
to move without colliding. Each particle is characterized by 




Fig. 1. Mechanical analog of a coupled oscillator network 
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a phase angle 6i G S 1 and has a preferred natural rotation 
frequency uji G R. Pairs of interacting particles i and j are 
coupled through an elastic spring with stiffness > 0. We 
refer to the Appendix [A] for a first principle modeling of the 
spring-interconnected particles depicted in Figure [T] 

Formally, each isolated particle is an oscillator with first- 
order dynamics 0$ = coj. The interaction among n such 
oscillators is modeled by a connected graph G(V,£,A) with 
nodes V = {l,...,n}, edges £ c V x V, and positive 
weights dij > for each undirected edge {i, j} G £ . Under 
these assumptions, the overall dynamics of the coupled 
oscillator network are 

0i = uji — 1 &ij sin(#i —Oj), i G {1, . . . , n} . (1) 

The rich dynamic behavior of the coupled oscillator model 
([T]) arises from a competition between each oscillator's 
tendency to align with its natural frequency Ui and the 
synchronization-enforcing coupling sin (6^ — 6j) with its 
neighbors. Intuitively, a weakly coupled and strongly het- 
erogeneous network does not display any coherent behavior, 
whereas a strongly coupled and sufficiently homogeneous 
network is amenable to synchronization, where all frequen- 
cies 0i(t) or even all phases 0i(t) become aligned. 

History, applications and related literature: The cou- 
pled oscillator model ([]]) has first been proposed by Arthur 
Winfree [4]. In the case of a complete interaction graph, 
the coupled oscillator dynamics ([]} are nowadays known as 
the Kuramoto model of coupled oscillators due to Yoshiki 
Kuramoto [5], [6]. Stephen Strogatz provides an excellent 
historical account in [7]. We also recommend the survey [8]. 

Despite its apparent simplicity, the coupled oscillator 
model ([]]) gives rise to rich dynamic behavior. This model 
is encountered in various scientific disciplines ranging from 
natural sciences over engineering applications to social net- 
works. The model and its variations appear in the study if 
biological synchronization phenomena such as pacemaker 
cells in the heart [9], circadian rhythms [10], neuroscience 
[1 1]— [13], metabolic synchrony in yeast cell populations 
[14], flashing fireflies [15], chirping crickets [16], biological 
locomotion [17], animal flocking behavior [18], fish schools 
[19], and rhythmic applause [20], among others. The coupled 
oscillator model ([T]) also appears in physics and chemistry in 
modeling and analysis of spin glass models [21], [22], flavor 
evolutions of neutrinos [23], coupled Josephson junctions 
[24], and in the analysis of chemical oscillations [25]. 

Some technological applications of the coupled oscillator 
model ([]} include deep brain stimulation [26], [27], vehicle 
coordination [19], [28]— [31], carrier synchronization without 



phase-locked loops [32], semiconductor lasers [33], [34], mi- 
crowave oscillators [35], clock synchronization in decentral- 
ized computing networks [36]-[41], decentralized maximum 
likelihood estimation [42], and droop-controlled inverters in 
microgrids [43]. Finally, the coupled oscillator model (p} 
also serves as the prototypical example for synchronization 
in complex networks [44] -[47] and its linearization is the 
well-known consensus protocol studied in networked control, 
see the surveys and monographs [48]-[50]. Various control 
scientists explored the coupled oscillator model |IJ as a 
nonlinear generalization of the consensus protocol [51]— [57]. 

Second-order variations of the coupled oscillator model 
(p} appear in synchronization phenomena, in population of 
flashing fireflies [58], in particle models mimicking animal 
flocking behavior [59], [60], in structure-preserving power 
system models, [61], [62] in network-reduced power system 
models [63], [64], in coupled metronomes [65], in pedes- 
trian crowd synchrony on London's Millennium bridge [66], 
and in Huy gen's pendulum coupled clocks [67]. Coupled 
oscillator networks with second-order dynamics have been 
theoretically analyzed in [8], [68]-[74], among others. 

Coupled oscillator models of the form ([T]) are also studied 
from a purely theoretic perspective in the physics, dynamical 
systems, and control communities. At the heart of the cou- 
pled oscillator dynamics is the transition from incoherence 
to synchrony. Here, different notions and degrees of synchro- 
nization can be distinguished [74]-[76], and the (apparently) 
incoherent state features rich and largely unexplored dynam- 
ics as well [47], [77]-[79]. In this article we will be particu- 
larly interested in phase and frequency synchronization when 
all phases Oi (t) become aligned, respectively all frequencies 
0i(t) become aligned. We refer to [7], [8], [19], [28], [31], 
[52], [53], [56], [64], [74]-[76], [80]-[95], [95]-[114] for an 
incomplete overview concerning numerous recent research 
activities. We will review some of literature throughout the 
paper and refer to the surveys [7], [8], [44]-[46], [74] for 
further applications and numerous additional theoretic results 
concerning the coupled oscillator model ([T]). 

Contributions and contents: In this paper, we introduce 
the reader to synchronization in networks of coupled oscil- 
lators. We present a sample of important analysis concepts 
in a tutorial style and from a control- theoretic perspective. 

In Section [TTJ we will review a set of selected technological 
applications which are directly tied to the coupled oscillator 
model (p} and also relevant to control systems. We will cover 
vehicle coordination, electric power networks, and clock 
synchronization in depth, and also justify the importance 
of the coupled oscillator model (p} as a canonical model. 
Prompted by these applications, we review the existing re- 
sults concerning phase synchronization, phase balancing, and 
frequency synchronization, and we also present some novel 
results on synchronization in sparsely-coupled networks. 



Section [IV] presents a collection of important results 
regarding phase synchronization, phase balancing, and fre- 
quency synchronization. By now the analysis methods for 
synchronization have reached a mature level, and we present 
simple and self-contained proofs using a sample of different 
analysis methods. In particular, we present one result on 
phase synchronization and one result on phase balancing in- 
cluding estimates on the exponential synchronization rate and 



the region of attraction (see Theorem |4.3| and Theorem |4.4| ). 
We also present some implicit and explicit, and necessary 
and sufficient conditions for frequency synchronization in 
the classic homogeneous case of a complete and uniformly- 
weighted coupling graphs (see Theorem |4.5| ). Concerning 
frequency synchronization in sparse graphs, we present two 
partially new synchronization conditions depending on the 



In particular, Section [III] introduces the reader to dif- 
ferent synchronization notions, performance metrics, and 
synchronization conditions. We illustrate these results with 
a simple yet rich example that nicely explains the basic 
phenomenology in coupled oscillator networks. 



algebraic connectivity (see Theorem |4.6| and Theorem |4.7| ). 

In our technical presentation, we try to strike a balance 
between mathematical precision and removing unnecessary 
technicalities. For this reason some proofs are reported in the 
appendix and others are only sketched here with references 
to the detailed proofs elsewhere. Hence, the main technical 
ideas are conveyed while the tutorial value is maintained. 

Finally, Section |V| concludes the paper. We summarize the 
limitations of existing analysis methods and suggest some 
important directions for future research. 

Preliminaries and notation: The remainder of this sec- 
tion introduces some notation and recalls some preliminaries. 

Vectors and functions: Let l n and n be the n-dimensional 
vector of unit and zero entries, and let 1^ be the orthogonal 
complement of l n in R n , that is, 1^ = {x G R n : x JL l n }. 
Given an n-tuple (xi, . . . , x n ), let x G R n be the associated 
vector with maximum and minimum elements x max and x m [ n . 
For an ordered index set 1 of cardinality \X\ and an one- 
dimensional array {xi}i e x, let diag({c$}iex) G Rl x l x l x l be 
the associated diagonal matrix. Finally, define the continuous 
function sine : R — » R by sinc(x) = sm(x)/x for 

Geometry on the n-torus: The set S 1 denotes the unit 
circle, an angle is a point G S 1 , and an arc is a connected 
subset of S 1 . The geodesic distance between two angles 
01, 62 G S 1 is the minimum of the counter-clockwise and 
the clockwise arc lengths connecting 0i and 02. With slight 
abuse of notation, let |0i — 02 1 denote the geodesic distance 
between two angles 0i,02 G S 1 . The n-torus is the product 
set T n = S 1 x • • • x S 1 is the direct sum of n unit circles. 
For 7 G [0, 2tt[, let Arc n (7) C T n be the closed set of angle 
arrays = (0i, . . . , n ) with the property that there exists 
an arc of length 7 containing all 0i, . . . , n . Thus, an angle 
array G Arc n (7) satisfies max^^j^ _ >n } \9i — 0j\ < 7. 
Finally, let Arc n (7) be the interior of the set Arc n (7). 

Algebraic graph theory: Let G(V, £ , A) be an undirected, 
connected, and weighted graph without self-loops. Let A G 
R nxn be its symmetric nonnegative adjacency matrix with 
zero diagonal, an — 0. For each node i G {1, . . . , n}, define 
the nodal degree by deg^ = J]j=i a ir Let L G R nxn be 
the Laplacian matrix defined by L = diag({deg i }^ =1 ) — A. 
If a number i G {1, . . . , |£|} and an arbitrary direction is 
assigned to each edge {i,j} G £, the (oriented) incidence 



matrix B G R nX l £: l is defined component- wise by Bm = 1 
if node k is the sink node of edge t and by Bm = — 1 if 
node k is the source node of edge t\ all other elements are 
zero. For x G M n , the vector B T x has components Xi — Xj 
corresponding to the oriented edge from j to i, that is, B T 
maps node variables Xi, Xj to incremental edge variables 
Xi — Xj. If diag({a^}^ 5j } G ^) is the diagonal matrix of edge 
weights, then L — B dmg{{aij}^^ e g)B T . If the graph 
is connected, then Ker (B T ) = Ker (L) = span(l n ), all 
n — 1 non-zero eigenvalues of L are strictly positive, and 
the second-smallest eigenvalue A2 (L) is called the algebraic 
connectivity and is a spectral connectivity measure. 

II. Applications of Kuramoto Oscillators 
Relevant to Control Systems 

The mechanical analog in Figure [T] provides an intuitive 
illustration of the coupled oscillator dynamics ([I]), and we 
reviewed a wide range of examples from physics, life sci- 
ences, and technology in Section |T| Here, we detail a set 
of selected technological applications which are relevant to 
control systems scientists. 

A. Flocking, Schooling, and Planar Vehicle Coordination 

An emerging research field in control is the coordination of 
autonomous vehicles based on locally available information 
and inspired by biological flocking phenomena. Consider a 
set of n particles in the plane R 2 , which we identify with 
the complex plane C. Each particle i G V = {l,...,n} 
is characterized by its position V{ G C, its heading angle 
Oi G S 1 , and a steering control law Ui(r,6) depending on 
the position and heading of itself and other vehicles. For 
simplicity, we assume that all particles have constant and 
unit speed. The particle kinematics are then given by [115] 

r 4 = e w * 1 

, ' ie{l,...,n}, (2) 

where i = is the imaginary unit. If the control ui is 

identically zero, then particle i travels in a straight line with 
orientation 0^(0), and if U{ = uji G R is a nonzero constant, 
then the particle traverses a circle with radius l/|o^|. 

The interaction among the particles is modeled by a 
possibly time- varying interaction graph G(V, £(£), A(t)) de- 
termined by communication and sensing patterns. Some 
interesting motion patterns emerge if the controllers use only 
relative phase information between neighboring particles, 
that is, Ui = uJo(t) + fi{6i — 6j) for {i, j} G £(t) and 
ujq : R> —> R. For example, the control ui = ouo(t) — 
K • Y^j=i a ij(t) sm (0i — 0j) with g am ^ G R results in 

En 
o^sin^-fl,), ieV. (3) 

The controlled phase dynamics ^ correspond to the coupled 
oscillator model ([TJ with a time-varying interaction graph 
with weights K • (t) and identically time- varying natural 
frequencies uji = c^o (t) for alH G {1, . . . , n}. The controlled 
phase dynamics ([3} give rise to very interesting coordination 
patterns that mimic animal flocking behavior [18] and fish 



schools [19]. Inspired by these biological phenomena, the 
controlled phase dynamics ^ and its variations have also 
been studied in the context of tracking and formation con- 
trollers in swarms of autonomous vehicles [19], [28]— [31]. 
A few trajectories are illustrated in Figure [2j and we refer to 
[19], [28]— [31] for other control laws and motion patterns. 

In the following sections, we will present various tools to 
analyze the motion patterns in Figure [2j which we will refer 
to as phase synchronization and phase balancing. 




Fig. 2. Illustration of the controlled dynamics j2j-{3j with n = 6 particles, 
a complete interaction graph, and identical and constant natural frequencies 
LJo(t) = 1, where K = 1 in panel (a) and K = — 1 in panel (b). The 
arrows depict the orientation, the dashed curves show the long-term position 
dynamics, and the solid curves show the initial transient position dynamics. 
It can be seen that even for this simple choice of controller, the resulting 
motion results in "synchronized" or "balanced" heading angles for K = ±1. 

B. Power Grids with Synchronous Generators and Inverters 

Here, we present the structure -pre serving power network 
model introduced in [61] and refer to [62, Chapter 7] for 
detailed derivation from a higher order first principle model. 
Additionally, we equip the power network model with a set 
of inverters and refer to [43] for a detailed modeling. 

Consider an alternating current (AC) power network mod- 
eled as an undirected, connected, and weighted graph with 
node set V = {1, . . . , n}, transmission lines £ C V x V, and 
admittance matrix Y = Y T G C nxn . For each node, consider 
the voltage phasor V* = |^|e l6> " corresponding to the phase 
Oi G S 1 and magnitude \Vi\ > of the sinusoidal solution 
to the circuit equations. If the network is lossless, then the 
active power flow from node i to j is sm(6i — 0j), where 
we used the shorthand a^- = \Vi\ • \Vj\ • $s(Yij). 

In the following, we assume that the node set is partitioned 
as V = Vi U V2 U V3, where V\ are load buses, V2 are con- 
ventional synchronous generators, and V3 are grid-connected 
direct current (DC) power sources, such as solar farms. The 
active power drawn by a load i e Vi consists of a constant 
term > and a frequency-dependent term Di$i with 
Di > 0. The resulting power balance equation is 

Di0i + I\ ti = -^2 aijSmiOi-Oj), zeVi. (4) 

If the generator reactances are absorbed into the admittance 
matrix, then the swing dynamics of generator i G V2 are 

En 
aijBhi(6i-6j), i G V 2 , (5) 



where Oi G S 1 and Oi G R 1 are the generator rotor angle 
and frequency, P m ^ > is the mechanical power input, and 
Mi > 0, and Di > are the inertia and damping coefficients. 

We assume that each DC source is connected to the AC 
grid via an DC/AC inverter, the inverter output impendances 
are absorbed into the admittance matrix, and each inverter is 
equipped with a conventional droop-controller. For a droop- 
controlled inverter i G V3 with droop- slope 1/Di > 0, the 
deviation of the power output Y^Jj=i a ij sm (^ — ®j) f rom 
its nominal value P<j,i > is proportional to the frequency 
deviation DiOi. This gives rise to the inverter dynamics 

En 
aijSmtfi-Oj), ieV 3 . (6) 

These power network devices are illustrated in Figure [3] 
Finally, we remark that different load models such as con- 
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Fig. 3. Illustration of the power network devices as circuit elements. Sub- 
figure (a) shows a transmission element connecting nodes i and j, Sub figure 
(b) shows a frequency-dependent load, Subfigure (c) shows an inverter con- 
trolled according to J5J, and Subfigure (d) shows a synchronous generator. 



stant power/current/susceptance loads and synchronous mo- 
tor loads can be modeled and analyzed by the same set of 
equations gJ-<5J, see [62]-[64], [116], [117]. 

Synchronization is pervasive in the operation of power 
networks. All generating units of an interconnected grid must 
remain in strict frequency synchronism while continuously 
following demand and rejecting disturbances. Notice that, 
with exception of the inertial terms Midi and the possibly 
non-unit coefficients Di, the power network dynamics (|4])-([6]) 
are a perfect electrical analog of the coupled oscillator model 
([T]) with uo = (—Pi,i : P m ,i,Pd,i)' Thus, it is not surprising 
that scientists from different disciplines recently advocated 
coupled oscillator approaches to analyze synchronization in 
power networks [43], [64], [69], [97], [114], [118]-[122]. 

The theoretic tools presented in the following sections 
establish how frequency synchronization in power networks 
depend on the nodal parameters P m ,i, Pd,i) as well 
as the interconnecting electrical network with weights a^. 
Ultimately, this deep understanding of synchrony gives us 
the correct intuition to design controllers and remedial action 
schemes preventing the loss of synchrony. 



C. Clock Synchronization in Decentralized Networks 

Another emerging technological application of the coupled 
oscillator model ([]} is clock synchronization in decentralized 
computing networks, such as wireless and distributed soft- 
ware networks. A natural approach to clock synchronization 
is to treat each clock as a coupled oscillator and follow a 
diffusion-based protocol to synchronize them, see the historic 
and recent surveys [36], [37], the landmark paper [38], and 
the interesting recent results [39]-[41]. 

Consider a set of distributed processors V = {1, . . . , n} 
interconnected in a (possibly directed) communication net- 
work. Each processor is equipped with an internal software 
clock, and these clocks need to be synchronized for dis- 
tributed computing and network routing tasks. For simplicity, 
we consider only analog clocks with continuous coupling 
since digital clocks are essentially discretized analog clocks 
and pulse-coupled clocks can be modeled continuously after 
a phase reduction and averaging analysis. 

For our purposes, the clock of processor i is a voltage- 
controlled oscillator which outputs a harmonic waveform 
Si(t) = sin(6^(t)), where 6i(t) is the accumulated instanta- 
neous phase. For uncoupled nodes, the phase 6i(t) evolves as 



i(0) 



2tt 



t I mod(27r) , i G {1, . . . , n} . 



where T nom > is the nominal period, Ti G R is an offset 
(frequency offset or skew), and 0$(O) G S 1 is the initial 
phase. To synchronize their internal clocks, the processors 
follow a diffusion-based protocol. In a first step, neighboring 
oscillators continuously communicate their respective wave- 
forms Si(t) to another. Second, through a phase detector each 
node measures a convex combination of phase differences as 

cvxi(0(t)) = =1 "ijWiW - OjV)) » * e {1, • • • , n} , 



where a^ > are convex (J^j 



1) and detector- 



specific weights, and / : S 1 — >> R is an odd 2 ^-periodic func- 
tion. Finally, cvx^ (#(£)) is fed to a (first-order and constant) 
phase-locked loop filter K whose output drives the local 
phase according to 



h{t) = 



2tt 

T). 



K -cvxMt)), *e{l,...,n}. (7) 



The goal of the synchronization protocol ([7]) is to synchronize 
the frequencies 6{ (t) or even the phases 6^ (t) in the processor 
network. For an undirected communication protocol, sym- 
metric weights = a^, and a sinusoidal coupling function 
/(•) = sin(-), the synchronization protocol ([7]) equals again 
the coupled oscillator model ([T]). 

The tools developed in the next section will enable us to 
state conditions when the protocol ^ successfully achieves 
phase or frequency synchronization. Of course, the protocol 
([7]) is merely a starting point, more sophisticated phase- 
locked loop filters can be constructed to enhance steady-state 
deviations from synchrony, and communication and phase 
noise as well as time-delays can be considered in the design. 



D. Canonical Coupled Oscillator Model 

The importance of the coupled oscillator model ([T]) does 
not stem only from the various examples listed in Sections [I] 
and|n| Even though model (p} appears to be quite specific (a 
phase oscillator with constant driving term and continuous, 
diffusive, and sinusoidal coupling), it is the canonical model 
of coupled limit-cycle oscillators [123]. In the following, 
we briefly sketch how such general models can be re- 
duced to model ([I]). We schematically follow the approaches 
[124, Chapter 10], [125] developed in the computational 
neuroscience community without aiming at mathematical 
precision, and we refer to [123], [126] for further details. 

Consider an oscillator modeled as a dynamical system with 
state x G M m and nonlinear dynamics x = f(x), which 
admit a locally exponentially stable periodic orbit 7 C W 71 
with period T > 0. By a change of variables, any trajectory 
in a local neighborhood of 7 can be characterized by a phase 
variable ip G S 1 with dynamics cp = Q, where ft = 2ir/T. 

Now consider a weakly forced oscillator of the form 



x = f(x) + e • S(t) , 



(8) 



where e > is sufficiently small and S(t) is a time-dependent 
forcing term. For small forcing eS(t), the attractive limit 
cycle 7 persists, and the phase dynamics are obtained as 

^ = n + eQ(^)5(t) + 0(6 2 ) J 

where Q((p) is the infinitesimal phase response curve (or 
linear response function), and we dropped higher order terms. 

Now consider n such limit cycle oscillators, where xi G 
R m is the state of oscillator i with limit cycle 7^ C lR m and 
period Ti > 0. We assume that the oscillators are weakly 
coupled with interaction graph G(V,£) and dynamics 

i i = /iW+ e Er- -i e 9ij(xi,Xj), i G {1, . . . , n} , (9) 

where gij(-) is the coupling function for the pair {i, j} G £. 
The coupling gij(-) can possibly be impulsive. The weak 
coupling in ^ can be identified with the weak forcing in 
([5]), and a transformation to phase coordinates yields 

(pi = fi t + eV .^ P Qi (%i (<Pi) , Xj (ipj)) , 

where fti — 2ir/Ti. The local change of variables 6i{t) — 
(Pi(t) — flit then yields the coupled phase dynamics 

0i = eJ2 f . QiiOi^ni^gijixiiOi^n^.XjiOj^njt)). 

An averaging analysis applied to the ^-dynamics results in 

0i = euji + e V hijiPi - Oj) , (10) 

where = /i^(0) and the averaged coupling functions are 



Mx) 



lim — 



i (Q i r)g ij (x i (ft i r),x j (ft j r-x))dT. 



Notice that the averaged coupling functions are 2tt- 
periodic and the coupling is diffusive. If all functions 
are odd, a first-order Fourier series expansion of yields 
hij(-) ~ a^ sin(-) as first harmonic with some coefficient 



a^. In this case, the dynamics fTO] ) in the slow time scale 
r = et reduce exactly to the coupled oscillator model (p}. 

This analysis justifies calling the coupled oscillator model 
([T]) the canonical model for coupled limit-cycle oscillators. 

III. Synchronization Notions and Metrics 

In this section, we introduce different notions of syn- 
chronization. Whereas the first four subsections address the 
commonly studied notions of synchronization associated 
with a coherent behavior and cohesive phases, Subsection 



III-E addresses the converse concept of phase balancing. 



A. Synchronization Notions 

The coupled oscillator model ([T]) evolves on T n , and 
features an important symmetry, namely the rotational in- 
variance of the angular variable 0. This symmetry gives rise 
to the rich synchronization dynamics. Different levels of syn- 
chronization can be distinguished, and the most commonly 
studied notions are phase and frequency synchronization. 

Phase synchronization: A solution 6 : M>o — » T n to the 
coupled oscillator model ([]} achieves phase synchronization 
if all phases 0i(t) become identical as t — » 00. 

Phase cohesiveness: As we will see later, phase synchro- 
nization can occur only if all natural frequencies uji are 
identical. If the natural frequencies are not identical, then 
each pairwise distance \0i(t) — 0j(t)\ can converge to a 
constant but not necessarily zero value. The concept of phase 
cohesiveness formalizes this possibility. For 7 G [0,7r[, let 
Ag(7) C T n be the closed set of angle arrays . . . , 6 n ) 
with the property \0i~9j \ < 7 for all {z, j} G £, that is, each 
pairwise phase distance is bounded by 7. Also, let Ag(7) be 
the interior of Ag(j). Notice that Arc n (7) C Ag(j) but the 
two sets are generally not equal. A solution 6 : R> — » T n 
is then said to be phase cohesive if there exists a length 
7 G [0,tt[ such that 6(t) G A g (t) for all t > 0. 

Frequency synchronization: A solution : M>o — » T n 
achieves frequency synchronization if all frequencies 0i{t) 



converge to a common frequency u, 
explicit synchronization frequency u. 



sync 



sync 



as t —> 00. The 
£ of the coupled 

oscillator model ([T]) can be obtained by summing over all 
equations in ([T]) as J27=i^ ~ E?=i ^* ^ n me frequency- 
synchronized case, this sum simplifies to J27=i ^sync = 
Y^i=i u i- ^ n conclusion, if a solution of the coupled oscillator 
model ([T]) achieves frequency synchronization, then it does so 
with synchronization frequency equal to u sync = Yl^i^i/ 71 - 
By transforming to a rotating frame with frequency co> sync and 
by replacing uoi by uji — co> sync , we obtain uj sync = (or equiv- 
alently w G In what follows, without loss of generality, 
we will sometimes assume that lj G 1^ so that lj, 



sync 



0. 



Remark 1 (Terminology): Alternative terminologies for 
phase synchronization include full, exact, or perfect synchro- 
nization. For a frequency- synchronized solution all phase 
distances \0i(i) — 0j(t)\ are constant in a rotating coordinate 
frame with frequency cj sync , and the terminology phase 
locking is sometimes used instead of frequency synchroniza- 
tion. Other commonly used terms include frequency locking, 
frequency entrainment, or also partial synchronization. □ 



Synchronization: The main object under study in most 
applications and theoretic analyses are phase cohesive and 
frequency- synchronized solutions, that is, all oscillators ro- 
tate with the same synchronization frequency, and all their 
pairwise phase distances are bounded. In the following, we 
restrict our attention to synchronized solutions with suffi- 
ciently small phase distances \6i—6j\ < 7 < tt/2 for {i, j} G 
£. Of course, there may exist other possible solutions, but 
these are not necessarily stable (see our analysis in Section 
[Tv| ) or not relevant in most application^] We say that a 
solution : R>o —> T n to the coupled oscillator model 
([T]) is synchronized if there exists sync G Ac (7) for some 
7 G [0,7r/2[ and uj sync G R (identically zero for uj G 1^) 



such that 0(t) = ; 



sync 



■ UJ. 



sync 



l n t (mod 2tt) for all t > 0. 



Synchronization manifold: The geometric object under 
study in synchronization is the synchronization manifold. 
Given a point r G S 1 and an angle s G [0, 2tt], let rot s (r) G 
S 1 be the rotation of r counterclockwise by the angle s. For 
(7*1, . . . , r n ) G T n , define the equivalence class 

[(n, . . . , r n )} = {(rot a (n), . . . , rot a (r n )) G T n | s G [0, 2tt]}. 

Clearly, if (ri,...,r n ) G Ag(7) for some 7 G [0,7r/2[, 
then [(n, . . . , r n )] C Ag(7). Given a synchronized solution 
characterized by sync G Ag(7) for some 7 G [0,7r/2[, 
the set [# sy nc] C Ac (7) is a synchronization manifold of 
the coupled-oscillator model (p}. Note that a synchronized 
solution takes value in a synchronization manifold due to 
rotational symmetry, and for uj G 1^- (implying cj sync = 0) a 
synchronization manifold is also an equilibrium manifold of 
the coupled oscillator model ([T]). These geometric concepts 
are illustrated in Figure [4] for the two-dimensional case. 



[6* 




A g (tt/2) 



Fig. 4. Illustration of the state space T 2 , the set Ag(tt/2), the synchro- 
nization manifold [6*] associated to a phase-synchronized angle array 0* = 
, #2) G Ag(0), and the tangent space with translation vector I2 at 0*. 



B. A Simple yet Illustrative Example 

The following example illustrates the different notions 
of synchronization introduced above and points out various 
important geometric subtleties occurring on the compact state 
space T 2 . Consider n = 2 oscillators with UJ2 > > uj\ = 
—uj2 • We restrict our attention to angles contained in an open 
half-circle: for angles 0i, 62 with |02 — 0i| < tt, the angular 

l For example, in power network applications the coupling terms 
dij sin(0i — Qj) are power flows along transmission lines {i,j} G £, and 
the phase distances \0i — O A are bounded well below tt/2 due to thermal 
constraints. In Subsection |III-E| we present a converse synchronization 
notion, where the goal is to maximize phase distances. 



difference 62 — 0i is the number in ] — 7r, tt[ with magnitude 
equal to the geodesic distance |#2 — #i| and with positive 
sign if and only if the counter-clockwise path length from 

01 to 62 on T 1 is smaller than the clockwise path length. 
With this definition the two-dimensional oscillator dynamics 
(01,62) can be reduced to the scalar difference dynamics 

02 — 01- After scaling time as t \-> t(uj2—uJi) and introducing 
k = 2ai2/(uJ2 — wi) the difference dynamics are 



d 

di 



(02 - 0i) = U0 2 - 0i ) := 1 - ^sin(0 2 - 0i) . (11) 



The scalar dynamics ( [TT] ) can be analyzed graphically by 
plotting the vector field /^(02 — 0i) over the difference 
variable 02 — 0i, as in Figure [5(a) Figure 5(a) displays a 
saddle-node bifurcation at k = 1. For k < 1 no equilibrium 
of ( [TT] ) exists, and for k > 1 an asymptotically stable 
equilibrium s tabie = arcsin(^ _1 ) G ]0,7r/2[ together with 
a saddle point sa ddie = arcsin(^ _1 ) G ]7r/2,7r[ exists. 

For 0(0) G Arc n (|0 sa ddie|) all trajectories converge ex- 
ponentially to s tabie, that is, the oscillators synchronize 
exponentially. Additionally, the oscillators are phase cohesive 
if an only if 0(0) G Arc n (|0 sa ddie|), where all trajectories 
remain bounded. For 0(0) Arc n (|0 sa ddi e |) the difference 
02 (t) — 0i(t) will increase beyond tt, and by definition 
will change its sign since the oscillators change orientation. 
Ultimately, 02 (t) — 0\ (t) converges to the equilibrium sta bie 
in the branch where 02 — 0i < 0. In the configuration space 
T 2 this implies that the distance |02^) — 0i(t)\ increases to 
its maximum value 7r and shrinks again, that is, the oscil- 
lators are not phase cohesive and revolve once around the 
circle before converging to the equilibrium manifold. Since 
sin(0 stab i e ) = sin(0 sad die) = strongly coupled oscillators 
with k ^> 1 practically achieve phase synchronization from 
every initial condition in an open semi-circle. In the critical 
case, k = 1, the saddle equilibrium manifold at tt/2 is 
globally attractive but not stable. An representative trajectory 



is illustrated in Figure 5(b) 





(a) Vector field {TTJ for 2 -0i > (b) Trajectory 0(t) for k = 1 

Fig. 5. Plot of the vector field {TTJ for various values of k, and a trajectory 
6{t) G T 2 for the critical case k = 1, where the dashed line is the saddle 
equilibrium manifold and ■ and • depict 0(0) and lim^oo 0(t). The non- 
smoothness of the vector field /(02 — 0\) at the boundaries {0,7r} is an 
artifact of the non-smoothness of the geodesic distance on T 2 



In conclusion, the simple but already rich 2 -dimensional 
case shows that two oscillators are phase cohesive and 
synchronize if and only if k > 1, that is, if and only if the 
coupling dominates the non-uniformity as 2ai2 > UJ2 — w\. 
The ratio 1/n determines the ultimate phase cohesiveness as 



well as the set of admissible initial conditions. For k ^> 1, 
practical phase synchronization is achieved for all angles 
in an open semi-circle. More general coupled oscillator 
networks display the same phenomenology, but the threshold 
from incoherence to synchrony is generally unknown. 

C. Synchronization Metrics 

The notion of phase cohesiveness can be understood as a 
performance measure for synchronization and phase synchro- 
nization is simply the extreme case of phase cohesiveness 
with lim^oo 6(t) G A G (0) = Arc n (0). An alternative 
performance measure is the magnitude of the so-called order 
parameter introduced by Kuramoto [5], [6]: 

The order parameter is the centroid of all oscillators repre- 
sented as points on the unit circle in C 1 . The magnitude r 
of the order parameter is a synchronization measure: if all 
oscillators are phase-synchronized, then r = 1, and if all os- 
cillators are spaced equally on the unit circle, then r = 0. The 
latter case is characterized in Subsection IIII-EI For a com- 
plete graph, the magnitude r of the order parameter serves as 
an average performance index for synchronization, and phase 
cohesiveness can be understood as a worst-case performance 
index. Extensions of the order parameter tailored to non- 
complete graphs have been proposed in [19], [52], [56]. 

For a complete graph and for 7 sufficiently small, the set 
A 0(7) reduces to Arc n (7), the arc of length 7 containing 
all oscillators. The order parameter is contained within 
the convex hull of this arc since it is the centroid of all 
oscillators, see Figure [6] In this case, the magnitude r of the 
order parameter can be related to the arc length 7. 




Fig. 6. Schematic illustration of an arc of length 7 £ [0, 7r], its convex hull 
(shaded), and the value • of the corresponding order parameter re 1 ^ with 
minimum magnitude r m [ n = cos (7/ 2) and maximum magnitude r max = 1. 

Lemma 3.1: (Shortest arc length and order parameter, 
[74, Lemma 2.1]) Given an angle array = (0i, . . . , n ) G 
T n with n > 2, let r(<9) = ±| YTj=i eWj I be the magnitude 
of the order parameter, and let 7(0) be the length of the 
shortest arc containing all angles, that is, G Arc n (7(0)). 
The following statements hold: 

1) if 7(g) G [0,tt], then r(0) G [cos( 7 (0)/2), 1]; and 

2) if G Arc n (7r), then 7(0) G [2 arccos(r(0)), tt]. 

D. Synchronization Conditions 

The coupled oscillator dynamics ([T]) feature (i) the syn- 
chronizing coupling described by the graph G(V,S,A) and 
(ii) the de- synchronizing effect of the non-uniform natural 



frequencies oj. Loosely speaking, synchronization occurs 
when the coupling dominates the non-uniformity. Various 
conditions have been proposed in the synchronization and 
power systems literature to quantify this trade-off. 

The coupling is typically quantified by the algebraic 
connectivity A 2 (L) [44], [45], [52], [64], [127], [128] or the 
weighted nodal degree deg^ = Xlj=i a u [64], [97], [117], 
[129], [130], and the non-uniformity is quantified by either 
absolute norms \\uj\\ p or incremental norms ||£> t cj|| p , where 
typically p G {2, 00}. Sometimes, these conditions can be 
evaluated only numerically since they are state-dependent 
[127], [129] or arise from a non-trivial linearization process, 
such as the Master stability function formalism [44], [45], 
[131]. In general, concise and accurate results are known only 
for specific topologies such as complete graphs [74], linear 
chains [108], and bipartite graphs [82] with uniform weights. 

For arbitrary coupling topologies only sufficient conditions 
are known [52], [64], [127], [129] as well as numerical 
investigations for random networks [89], [98], [99], [128], 
[132]. Simulation studies indicate that these conditions are 
conservative estimates on the threshold from incoherence to 
synchrony. Literally, every review article on synchronization 
draws attention to the problem of finding sharp synchroniza- 
tion conditions [7], [8], [44]-[46], [74], [114]. 

E. Phase Balancing and Splay State 

In certain applications in neuroscience [1 1]— [13], deep- 
brain stimulation [26], [27], and vehicle coordination [19], 
[28]— [31], one is not interested in the coherent behavior 
with synchronized (or nearly synchronized) phases, but rather 
in the phenomenon of synchronized frequencies and de- 
sychronized phases. 

Whereas the phase- synchronized state is characterized by 
the order parameter r achieving its maximal (unit) magni- 
tude, we say that a solution : R> — » T n to the coupled 
oscillator model ([]} achieves phase balancing if all phases 
0i(t) converge to Bal n = {0 G T n : r(0) = | £ YTj=i eWj I = 
0} as t — > 00, that is, the oscillators are distributed over the 
unit circle S 1 , such that their centroid re 1 ^ vanishes. We refer 
to [28] for a geometric characterization of the balanced state. 

One balanced state of particular interest in neuroscience 
applications [1 1]— [13], [26], [27] is the so-called splay state 
{0 G T n : 6i = i-27r/n + (p (mod 2tt) , <p G 8 1 , i G 
{l,...,n}} C Bal n corresponding to phases uniformly 
distributed around the unit circle S 1 with distances 2ii/n. 
Other highly symmetric balanced states consist of multiple 
clusters of collocated phases, where the clusters themselves 
are arranged in splay state, see [28], [29]. 

IV. Analysis of Synchronization 

In this section we present several analysis approaches to 
synchronization in the coupled oscillator model ([T]). We begin 
with a few basic ideas to provide important intuition as well 
as the analytic basis for further analysis. 



A. Some Simple Yet Important Insights 

The potential energy U : T n — » R of the elastic spring 
network in Figure [T] is, up to an additive constant, given by 

u w = E , i}ef a ^ ( x - cos (^ - i)) ■ w 

By means of the potential energy, the coupled oscillator 
model ([]} can reformulated as the forced gradient system 



0i =Ui-ViU(0), i e {!,..., n} ; 



(13) 



where Vi~U{0) = ^-U(0) denotes the partial derivative. 
It can be easily verified that the phase- synchronized state 
Q i = 0- for all {i, j} G £ is a local minimum of the potential 
energy ( fT2| ). The gradient formulation ( [T3] ) clearly empha- 
sizes the competition between the synchronization-enforcing 
coupling through the potential U(0) and the synchronization- 
inhibiting heterogeneous natural frequencies uji. 

We next note that uj has to satisfy certain bounds, relative 
to the weighted nodal degree, in order for a synchronized 
solution to exist. 

Lemma 4.1: (Necessary sync conditions) Consider the 
coupled oscillator model ([T]) with graph G(V,£,A), fre- 
quencies uj G 1^, and nodal degree deg^ = Y^j=i a ij f° r 
each oscillator i G {1, . . . , n}. If there exists a synchronized 
solution G Ac (7) for some 7 G [0, 7r/2], then the 
following conditions hold: 

1) Absolute bound: For each node i G {1, . . . , n}, 



deg f sin(7) > \uJi 



(14) 



2) Incremental bound: For all distinct z, j G {1, . . . , n}, 

(deg^ + deg^-) sin(7) > \uJi - ujj | . (15) 

Proof: Since G 1^, the synchronization frequency 
is zero, and phase and frequency synchronized solutions 



^sync 



are equilibrium solutions determined by the equations 



G{l,...,n}. (16) 



Since sin(0» -0j) G [- sin(7), + sin(7)] for G Ac (7), the 
equilibrium equations ( fT6| ) have no solution if condition fl4| ) 
is not satisfied. Since uj G 1^, an incremental bound on uj 
seems to be more appropriate than an absolute bound. The 
subtraction of the ith and jth equation ([16]) yields 

En 
fc=l (a ik sin(6> i - k ) - a jk sin(0j - fc )) . 

Again, since the coupling is bounded, the above equation has 
no solution in Ac (7) if condition ( p3] ) is not satisfied. ■ 

The following result is fundamental for various approaches 
to phase and frequency synchronization. To the best of the 
authors' knowledge this result has been first established in 
[133], and it has been reproved numerous times. 

Lemma 4.2: (Stable synchronization in Ag(tt/2)) Con- 
sider the coupled oscillator model ([]} with a connected 
graph G(V,S,A) and frequencies uj G 1^. The following 
statements hold: 



1) Jacobian: The Jacobian J(0) of the coupled oscillator 
model ([]} evaluated at G T n is given by 

J{9) = -Bdmgdciij cos^ - j )} {iJ}ee )B T ; 

2) Local stability and uniqueness: If there exists an 
equilibrium 6* G Ag(tt/2), then 

(i) —J(Q*) is a Laplacian matrix; 

(ii) the equilibrium manifold [0*] G A g (tt/2) is 
locally exponentially stable; and 

(iii) this equilibrium manifold is unique in Ag(tt/2). 

Y^ =1 aijsm(0i - Oj)) = 



Proof: Since -J^ (uji 



- E™=i a ij cos (°i - Oj) and m~\ u i - T,™=i a ij s ' m (°i - 

Oj)) = cos(0i — Oj), we obtain that the Jacobian is 
equal to minus the Laplacian matrix of the connected graph 
G(V,£,A) with the (possibly negative) weights = 
cos(0i — Oj). Equivalently, in compact notation J(0) = 
— B diag({a^- cos(#i — Oj)}^jy e g)B T . This completes the 
proof of statement 1). 

The Jacobian J{0) evaluated for an equilibrium 0* G 
Ag(tt/2) is minus the Laplacian matrix of the graph 
G(V,£,A) with strictly positive weights a^ = a^ cos(0* — 
Oj) > for every {i,j} G £. Hence, J(0*) is negative 
semidefinite with the nullspace l n arising from the rotational 
symmetry, see Figure [4] Consequently, the equilibrium point 
0* G Ag(tt/2) is locally (transversally) exponentially sta- 
ble, or equivalently, the corresponding equilibrium manifold 
[0*] G Ag(7t/2) is locally exponentially stable. 

The uniqueness statement follows since the right-hand side 
of the coupled oscillator model (p} is a one-to-one function 
(modulo rotational symmetry) for G Ag(tt/2), see [134, 
Corollary 1]. This completes the proof of statement 2). ■ 



By Lemma 4.2 any equilibrium in Ac(7r/2) is stable 
which supports the notion of phase cohesiveness as a per- 
formance metric. Since the Jacobian J(0) is the negative 
Hessian of the potential U(0) defined in ( [T2| ), Lemma 4.2 
also implies that any equilibrium in Ag*(7t/2) is a local 
minimizer of U(0). Of particular interest are so-called S 1 - 
synchronizing graphs for which all critical points of ( [T2| ) 
are hyperbolic, the phase- synchronized state is the global 
minimum of U(0), and all other critical points are local 
maxima or saddle points. The class of S 1 -synchronizing 
graphs includes, among others, complete graphs and acyclic 
graphs [100]-[103]. 

These basic insights motivated various characterizations 
and explorations of the critical points and the curvature of 
the potential U(0) in the literature on synchronization [52], 
[64], [74], [89], [93], [100], [100]-[103], [103] as well as 
on power systems [61], [116], [127], [129], [133]-[137]. 

B. Phase Synchronization 

If all natural frequencies are identical, uj^ = uj for all i G 
{1, . . . , n}, then a transformation of the coupled oscillator 
model ([!} to a rotating frame with frequency uj leads to 

En 
^aijsm(0i -Oj), zG {!,..., n}. (17) 



The analysis of the coupled oscillator model ^Tf\ is par- 
ticularly simple and local phase synchronization can be 
concluded by various analysis methods. A sample of dif- 
ferent analysis schemes (by far not complete) includes 
the contraction property [54], [64], [92], [100], [138], 
quadratic Lyapunov functions [52], [64], linearization [81], 
[103], or order parameter and potential function arguments 
[28], [56], [80]. 

The following theorem on phase synchronization summa- 
rizes a collection of results originally presented in [28], [54], 
[56], [74], [100], [103], and it can be easily proved given the 



insights developed in Subsection IV-A 

Theorem 4.3: (Phase synchronization) Consider the cou- 
pled oscillator model ([]]) with a connected graph G(V,£,A) 
and with frequency uj G W 1 (not necessarily zero mean). The 
following statements are equivalent: 

(i) Stable phase sync: there exists a locally exponentially 
stable phase-synchronized solution 9 G Arc n (0) (or a 
synchronization manifold [9] G Ag(0)); and 

(ii) Uniformity: there exists a constant uj G R such that 
uji = uj for all i G {1, . . . , n}. 

If the two equivalent cases (i) and (ii) are true, the following 
statements hold: 

1) Global convergence: For all initial angles 9(0) G T n 
all frequencies 9i(t) converge to uj and all phases 
9i(t) — uot (mod 2tt) converge to the critical points 
{9 G T n : VU(9) = n }; 

2) Semi-global stability: The region of attraction of the 
phase-synchronized solution 9 G Arc n (0) contains the 
open semi-circle Arc n (7r), and each arc Arc n (7) is 
positively invariant for every arc length 7 < 7r; 

3) Explicit phase: For initial angles in an open semi- 
circle 9(0) G Arc n (7r), the asymptotic synchronization 
phase is given b)Q«W = Er=i^(0)/n+a;t (mod 2tt); 

4) Convergence rate: For every initial angle 9(0) G 
Arc n (7) with 7 < tt, the exponential convergence 
rate to phase synchronization is no worse than A ps = 
— \2(L) sinc(7); and 

5) Almost global stability: If the graph G(V,£,A) is 
S 1 -synchronizing, the region of attraction of the phase- 
synchronized solution 9 G Arc n (0) is almost all of T n . 



Proof: Implication (i) 



(ii): By assumption, there 



exist constants 



' sync 



and uj, 



sync 



such that 9i(t) 



#sync + ^sync^ (mod 2tt). In the phase- synchronized case, the 
dynamics ([]} then read as cj sync = uji for all i G {1, . . . , n}. 
Hence, a necessary condition for the existence of phase 
synchronization is that all uji are identical. 

Implication (ii) (i): Consider the model ([]} written 

in a rotating frame with frequency uj as in ^Tf\ . Note that the 
set of phase-synchron ized solutions Aq(0) is an equilibrium 
manifold. By Lemma 4.2 we conclude that Aq(0) is locally 
exponentially stable. This concludes the proof of (i) <^> (ii). 



2 This "average" of angles (points on S 1 ) is well-defined in an open 
semi-circle. If the parametrization of 6 has no discontinuity inside the arc 
containing all angles, then the average can be obtained by the usual formula. 



Statement 1 ): Note that ( fT7] ) can be written as the gradient 
flow 9 = —VU(9), and the corresponding potential function 
U(9) is non-increasing along trajectories. Since the sublevel 
sets of U(9) are compact and the vector field VU(9) is 
smooth, the invariance principle [139, Theorem 4.4] asserts 
that every solution converges to set of equilibria of §FJ\ . 

Statements 2): The coupled oscillator model $FI\ can be 
re-written as the consensus-type system 



j j=1 



b ij (9)-(9 i -9 j ), ie {!,..., n}, (18) 



where the weights bij(9) = a^ sinc(^ — 9j) depend ex- 
plicitly on the system state. Notice that for 9 G Arc n (7) 
and 7 < 7r the weights b^ (9) are upper and lower bounded 
as bij(9) G [aij sine (7), a^] Assume that the initial angles 
9i(0) belong to the set Arc n (7), that is, they are all contained 
in an arc of length 7 G [0, 7r[. In this case, a natural Lyapunov 
function to establish phase synchronization can be obtained 
from the contraction property, which aims at showing that 
the convex hull containing all oscillators is decreasing, see 
[54], [64], [92], [100], [140] and the review [138, Section 2]. 

Recall the geodesic distance between two angles on S 1 
and define the continuous function V : T n — >> [0, tt] by 



V(ip) =max{|^-^| I i,j G {!,..., n}}. 



(19) 



Notice that, if all angles are contained in an arc at time t, 
then the arc length V(9(t)) = max^ jG { l5 _ ?n } \9i(t) — 9j(t)\ 
is a Lyapunov function candidate for phase synchronization. 
Indeed, it can be shown that V(9(t)) decreases along trajec- 
tories of ( [T8] ) for 9(0) G Arc n (7) and for all 7 < tt. The 
analysis is complicated by the following fact: the function 
V(9(t)) is continuous but not necessarily differentiable when 
the maximum geodesic distance (that is, the right-hand side 
of fl9])), is attained by more than one pair of oscillators. We 
omit the explicit calculations here and refer to [54], [64], 
[74], [83], [92] for a detailed analysis. 

Statement 3): By statement 2), the set Arc n (7r) is pos- 
itively invariant, and for 9(0) G Arc n (7r) the average 
^i=i9i(t)/n is well defined for t > 0. A summation 
over all equations of the model $FJ\ yields J^Li^W = 
0, or equivalently, J^ILi^W is constant for all t > 0. 
In particular, for t = we have that ^ILi^W = 
J27=i ^i(O) anc * f° r a phase- synchronized solution we have 
that J27=i ^sync = J27=i ^(0)- Hence, the explicit synchro- 
nization phase is given by 5^=1 ^(0)/ n - ^ n me ori g ma l 
coordinates (non-rotating frame) the synchronization phase 
is given by J27=i ^(0)/ n + u ^ ( m °d 2tt). 

Statement 4): Given the invariance of the set Arc n (7) for 
any 7 < tt, the system ( fl"8] ) can be analyzed as a linear 
time- varying consensus system with initial condition 9(0) G 
Arc n (7), and bounded time- varying weights bij(9(t)) G 
[aij sine (7), a^] for all t > 0. The worst-case convergence 
rate A ps can then be obtained by a standard symmetric 
consensus analysis, see [52], [53], [64], [74]. For instance, 
it can be shown that the deviation of the angles 9(t) from 
their average, \\9(t) — (J27=i Q(t)/ n )ln\\2 ( me disagreement 
function) decays exponentially with rate A ps . 



Statement 5): By statement 1), all solutions of system flT] ) 
converge to the set of equilibria, which equals the set of 
critical points of the potential U (6). By the definition of S 1 - 
synchronizing graphs, the phase- synchronized equilibrium 
manifold Arc n (0) is the only stable equilibrium set, and all 
others are unstable. Hence, for all initial condition 0(0) G 
T n , which are not on the stable manifolds of unstable 
equilibria, the corresponding solution 6(t) will reach the 
phase-synchronized equilibrium manifold Arc n (0). ■ 

Remark 2: (Control-theoretic perspective on synchro- 



nization) As established in Theorem 4.3 the set of phase- 
synchronized solutions Arc n (0) of the coupled oscillator 
model (p} is locally stable provided that all natural frequen- 
cies are identical. For non-uniform (but sufficiently identical) 
natural frequencies, phase synchronization is not possible but 
a certain degree of phase cohesiveness can still be achieved. 
Hence, the coupled oscillator model ([T]) can be regarded as 
an exponentially stable system subject to the disturbance 
uj G 1^, and classic control-theoretic concepts such as input- 
to-state stability, practical stability, and ultimate boundedness 
[139] or their incremental versions [141] can be used to study 
synchronization. In control-theoretic terminology, synchro- 
nization and phase cohesiveness can then also be described as 
"practical phase synchronization". Compared to prototypical 
nonlinear control examples, various additional challenges 
arise in the analysis of the coupled oscillator model (p} due to 
the bounded and non-monotone sinusoidal coupling and the 
compact state space T n containing numerous equilibria; see 



the analysis approaches in Section IV and [64], [74], [95]. □ 



C. Phase Balancing 

In general, only few results are known about the phase 
balancing problem. This asymmetry is partially caused by 
the fact that phase synchrony is required in more applications 
than phase balancing. Moreover, the phase-synchronized set 
Arc n (0) admits a very simple geometric characterization, 
whereas the phase-balanced set Bal n has a complicated struc- 
ture consisting of numerous disjoint subsets. The number 
of these subsets grows with the number of nodes n in a 
combinatorial fashion. 

Consider the coupled oscillator model flT] ) with identical 
natural frequencies. By inverting the direction of time, we get 



En 
_ i a,ij sm(6i 



G {!,..., n}. (20) 



In the following, we say that the interaction graph G(V,£,A) 
is circulant if the adjacency matrix A = A T is a circulant 
matrix. Circulant graphs are highly symmetric graphs includ- 
ing complete graphs, bipartite graphs, and ring graphs]^] For 
circulant and uniformly weighted graphs, the coupled oscil- 
lator model pO} achieves phase balancing. The following 
theorem summarizes different results, which were originally 
presented in [28], [29]. 

Theorem 4.4: (Phase balancing) Consider the coupled 
oscillator model ([20]) with a connected, uniformly weighted, 



and circulant graph G(V,£,A). The following statements 
hold: 

1) Local phase balancing: The phase-balanced set Bal n 
is locally asymptotically stable; and 

2) Almost global stability: If the graph G(V,£,A) is 
complete, then the region of attraction of the stable 
phase-balanced set Bal n is almost all of T n . 

The proof of Theorem |4.4| follows a similar reasoning 
as the proof of Theorem |4.3| convergence is established 
by potential function arguments and local (in) stability of 
equilibria by Jacobian arguments. We omit the proof here 
and refer to [28, Theorem 1] and [29, Theorem 2] for details. 

For general connected graphs, the conclusions of Theorem 
4.4 are not true. As a remedy to achieve locally stable and 
globally attractive phase balancing, higher order models need 
to be considered, see the models proposed in [29], [56]. 

D. Synchronization in Complete Networks 

For a complete coupling graph with uniform weights = 
K/n, where K > is the coupling gain, the coupled oscil- 
lator model ([T]) reduces to the celebrated Kuramoto model 

^=^-^^J =i sin(^-^), ie{l,...,n}. (21) 

By means of the order parameter re 1 ^ = ^ Y^=i el9j > me 
Kuramoto model ([21} can be rewritten in the insightful form 



% = uji - Kr sin(0i - iji) , ie {!,..., n}. 



(22) 



Equation ( [22] ) gives the intuition that the oscillators syn- 
chronize by coupling to a mean field represented by the 
order parameter re 1 ^. Intuitively, for small coupling strength 
K each oscillator rotates with its natural frequency u^, 
whereas for large coupling strength K all angles Oi(t) will 
be entrained by the mean field re 1 ^ and the oscillators 
synchronize. The threshold from incoherence to synchroniza- 
tion occurs for some critical coupling identical- This phase 
transition has been the source of numerous investigations 
starting with Kuramoto's analysis [5], [6]. Various necessary, 
sufficient, implicit, and explicit estimates of the critical 
coupling strength identical for both the on-set as well as 
the ultimate stage of synchronization have been proposed 
[5]-[8], [28], [52], [53], [64], [74], [75], [82]-[87], [95]- 
[97], [100], [103]-[107], [110], and we refer to [74] for a 
comprehensive overview. 

The mean field approach to the equations ( [22] ) can be made 
mathematically rigorous by a time-scale separation [96] or 
in the continuum limit as the number of oscillators tends 
to infinity and the natural frequencies uj are sampled from 
a distribution function g : R — >> M>o. In the continuum 
limit and for a symmetric, continuous, and unimodal dis- 
tribution g(uj), Kuramoto himself showed in an insightful 
and ingenuous analysis [5], [6] that the incoherent state (a 
uniform distribution of the oscillators on the unit circle S 1 ) 
supercritically bifurcates for the critical coupling strength 



3 Further info on circulant graphs and a gallery can be found at 

|http : //mathworld . wolfram. com/ Circulant G raph . html 



-^critical — 



ng(0) 



(23) 



In [8], [87], [104], it was found that the bipolar (bimodal 
double-delta) distribution (respectively the uniform distri- 
bution) yield the largest (respectively smallest) threshold 
^critical over all distributions g(uu) with bounded support. We 
refer [7], [8] for further references and to [88], [109]-[111] 
for recent contributions on the continuum limit model. 

In the finite-dimensional case, the necessary synchroniza- 
tion condition ( [T5] ) gives a lower bound for identical as 

Ti 

K > — — • (o; max - uj min ) . (24) 

2[n — 1) 

Three recent articles [84]-[86] independently derived a 
set of implicit consistency equations for the exact critical 
coupling strength identical for which synchronized solutions 
exist. Verwoerd and Mason provided the following implicit 
formulae to compute identical [86, Theorem 3]: 

En i 

, (25) 

2]T yi-(o,/ u *)2 = ^. i/v/i-^/^)^ 

z — 'i=l z — 'i=l 

where ft* = ^ - cj sync and G [H^H^ , 2 H^H^]. The im- 
plicit formulae ([25]) can also be extended to bipartite graphs 
[82]. A local stability analysis is carried out in [84], [85]. 

From the point of analyzing or designing a sufficiently 
strong coupling, the exact formulae ( [25] ) have three draw- 
backs. First, they are implicit and thus not suited for perfor- 
mance or robustness estimates in case of additional coupling 
strength for a given K > identical- Second, the corresponding 
region of attraction of a synchronized solution is unknown. 
Third and finally, the particular natural frequencies uji are 
typically time-varying, uncertain, or even unknown in the 
applications listed in Section [I] In this case, the exact value 
of identical needs to be estimated in continuous time, or a 
conservatively strong coupling K ^> identical has to be chosen. 

The following theorem states an explicit bound on the crit- 
ical coupling strength together with performance estimates, 
convergence rates, and a guaranteed semi-global region of 
attraction for synchronization. This bound is tight and thus 
necessary and sufficient when considering arbitrary distribu- 
tions of the natural frequencies with compact support. The 
result has been originally presented in [74, Theorem 4.1]. 

Theorem 4.5: (Synchronization in the Kuramoto 
model) Consider the Kuramoto model ( [2T] ) with natural 
frequencies uo = (lji, . . . , oj n ) and coupling strength K. The 
following three statements are equivalent: 

(i) the coupling strength K is larger than the maximum 
non-uniformity among the natural frequencies, that is, 

K > -^critical ^max — ^min 5 

(26) 

(ii) there exists an arc length 7 max G ]7r/2, tt] such that the 
Kuramoto model pT| ) synchronizes exponentially for 
all possible distributions of the natural frequencies uoi 
supported on the compact interval [c^mm, cj max ] and for 
all initial phases 0(0) G Arc n (7 max ); and 

(iii) there exists an arc length 7 m i n G [0,7r/2[ such that the 
Kuramoto model ( [2T] ) has a locally exponentially stable 
synchronization manifold in Arc n (7 m i n ) for all possible 



distributions of the natural frequencies supported on 
the compact interval [cj m i n , o; max ]. 

If the three equivalent conditions (i), (ii), and (iii) hold, 
then the ratio identical and the arc lengths 7^ G [0, 7r/2[ 
and 7 max G ]7r/2,7r] are related uniquely via sin(7 m i n ) = 
sin(7 max ) = i^ cr i t i ca i/K, and the following statements hold: 

1) phase cohesiveness: the set Arc n (7) C Ag(7) is 
positively invariant for every 7 G [7min, 7max], and each 
trajectory starting in Arc n (7 max ) approaches asymptot- 
ically Arc n (7 min ); 

2) frequency synchronization: the asymptotic synchro- 
nization frequency is the average frequency cj sync = 
- Y^i=i ^i* anc *' gi yen phase cohesiveness in Arc n (7) 
for some fixed 7 < 7r/2, the exponential synchroniza- 
tion rate is no worse than Ak = — i\~cos(7); and 

3) order parameter: the asymptotic value of the mag- 
nitude of the order parameter, denoted by = 
lim^oo 1| YJj=i e Wj ^\, is bounded as 

1 > > cos(^p) = ^ + 

Proof: In the following, we sketch the proof of Theo- 
rem [43] and refer to [74, Theorem 4.1] for further details. 

Implication (i) =4> (ii): In a first step, it is shown that the 
phase cohesive set Arc n (7) is positively invariant for every 
7 £ [7min?7max]- By assumption, the angles 6i(t) belong to 
the set Arc n (7) at time t = 0, that is, they are all contained 
in an arc of length 7. We aim to show that all angles remain 
in Arc n (7) for all subsequent times t > by means of the 
contraction Lyapunov function ([19]). Note that Arc n (7) is 
positively invariant if and only if V(0(t)) does not increase at 
any time t such that V(6(t)) = 7. The upper Dini derivative 
of V(6(t)) along trajectories of ( [2T] ) is given by 

h^O h 

Written out in components and after trigonometric simplifi- 
cations [74], we obtain that the derivative is bounded as 

D+V{0{t)) < co max - cj min - Ksm(j) . 

It follows that the length of the arc formed by the angles is 
non-increasing in Arc n (7) if and only if 

Ksin(j) > identical, (27) 

where identical is as stated in equation ( [26] ). For 7 G [0, tt] the 
left-hand side of ( [27] ) is a concave function of 7 that achieves 
its maximum at 7* = tt/2. Therefore, there exists an open set 
of arc lengths 7 G [0, tt] satisfying equation ( [27] ) if and only 
if equation ( [27] ) is true with the strict equality sign at 7* = 
7r/2, which corresponds to condition ( [26] ). Additionally, if 
these two equivalent statements are true, then there exists 
a unique 7^ G [0,7r/2[ and a 7 max G ]7r/2,7r] that satisfy 
equation ( [27] ) with the equality sign, namely sin (7^) = 
sin(7 max ) = K CYitica i/K. For every 7 G [7min,7max] it follows 
that the arc-length V(0(t)) is non-increasing, and it is strictly 
decreasing for 7 G ]7min? 7max[- Among other things, this 



shows that statement (i) implies statement 1). By means of 
Lemma [XT] statement 3) then follows from statement 1). 

The frequency dynamics of the Kuramoto model ( [2T] ) can 
be obtained by differentiating the Kuramoto model ( |2T] ) as 



d 



(28) 



where a>ij(t) = (K/ri) cos(6i(t) — Oj(t)). For K > identical, 
we just proved that for every 6(0) G Arc n (7 max ) and for all 
7 G ]7mim7max[ there exists a finite time T > such that 
6(t) G Arc n (7) for all t >T. Consequently, the terms aij(t) 
are strictly positive for all t > T. Notice also that system 
( [28] ) evolves on the tangent space of T n , that is, the Euclidean 
space R n . Now fix 7 G ]7min,7r/2[ and let T > such that 
CLij (t) > for all t > T. In this case, the frequency dynamics 
( [28] ) can be analyzed as linear time-varying consensus sys- 
tem. Consider the disagreement vector x = 6 — cj sync l n as 
an error coordinate. By standard consensus arguments [48]- 
[50], it can be shown that the disagreement vector satisfies 
\\x(t)\\ < ||x(0)||e" AKt for alU > T. This proves statement 
2) and the implication (i) (ii). 

Implication (ii) => (i): To show that condition ( [26] ) is 
also necessary for synchronization, it suffices to construct a 
counter example for which K < ^critical and the oscillators 
do not achieve exponential synchronization even though all 
Ui G [cj min , cj max ] and (9(0) G Arc n (7) for every 7 G ]tt/2, tt]. 
A basic instability mechanism under which synchronization 
breaks down is caused by a bipolar distribution of the natural 
frequencies. Let the index set {1, . . . , n} be partitioned by 
the two non-empty sets X\ and X 2 . Let uoi = uo^ m for z Gli 
and 00 i = cj max for i G X2, and assume that at some time 
t > it holds that 0i(t) = -7/2 for i G X 1 and 0i(t) = 
+7/2 for i G X 2 and for some 7 G [0, 7r[. By construction, 
at time t all oscillators are contained in an arc of length 
7 G [0,7r[. Assume now that if < identical and the oscillators 
synchronize. It can be shown [74] that the evolution of the 
arc length V(0(t)) satisfies the equality 



D+V(8(t)) 



K sin (7) . 



(29) 



Clearly, for K < identical the arc length V(6(t)) = 7 is 
increasing for any arbitrary 7 G [0, tt]. Thus, the phases are 
not bounded in Arc n (7). This contradicts the assumption 
that the oscillators synchronize for K < identical from every 
initial condition 6(0) G Arc n (7). For K = if critical, we 
know from [84], [85] that phase-locked equilibria have a 
zero eigenvalue with a two-dimensional Jacobian block, and 
thus synchronization cannot occur. This instability via a 
two-dimensional Jordan block is also visible in ( [29] ) since 
D + V(6(t)) is increasing for 6(t) G Arc n (7), 7 G ]7r/2,7r] 
until all oscillators change orientation, just as in the example 
in Subsection |III-B [ This proves the implication (ii) (i). 



Equivalence (i),(ii) (Hi): The proof relies on Jacobian 
arguments and will be omitted here, see [74] for details. ■ 

Theorem |4.5| places a hard bound on the critical coupling 
strength identical for all distributions of uji supported on the 
compact interval [cj m i n , cj max ]. For a particular distribution 



g(uj) supported on [co> m i n , cj max ] the bound ( [26] ) is only suffi- 
cient and possibly a factor 2 larger than the necessary bound 
( |24| ). The exact critical coupling lies somewhere in between 
and can be obtained from the implicit equations ( [25] ). 

Since the bound ( [26] ) on -/^critical is exact [74] for the worst- 
case bipolar distribution uji G {o; m i n , cj max }, Figure |7] reports 
numerical findings for the other extreme case [87] of a uni- 
form distribution g(u) = 1/2 supported for uji G [—1,1]. All 




Fig. 7. Statistical analysis of the necessary and explicit bound {24) ( ), 
the exact and implicit bound {25} (O), and the sufficient, tight, and explicit 
bound {26) (□) for n G [2, 300] oscillators in a semi-log plot, where the 
coupling gains for each n are averaged over 1000 simulations. 

three displayed bounds are identical for n = 2 oscillators. As 
n increases, the sufficient bound ( [26] ) converges to the width 
^max — ^min = 2 of the support of g(uj), and the necessary 
bound ( [24] ) accordingly to half of the width. The exact bound 
( [25] ) converges to 4(cj max — cj m i n )/(27r) =4/7r in agreement 
with condition ( [23] ) predicted for the continuum limit. 

Finally, let us mention that Theorem |4.5| can be extended 
to discontinuously switching and slowly time- varying natural 
frequencies [74]. For a particular sampling distribution g(u), 
the critical quantity in condition ( [26] ), the support cj max — u m [ n , 
can be estimated by extreme value statistics, see [89]. 

E. Synchronization in Sparse Networks 

As summarized in Subsection |III-D| the quest for sharp 
and concise synchronization for non-complete coupling 
graph G(V,£,A) is an important and outstanding problem 
emphasized in every review article on coupled oscillator 
networks [7], [8], [44]-[46], [74]. The approaches known for 
phase synchronization in arbitrary graphs or the contraction 
approach to frequency synchronization (used in the proof of 
Theorem [43] ) do not generally extend to arbitrary natural fre- 
quencies wGl^ and connected coupling graphs G(V, £ , A), 
or do so only under extremely conservative conditions. 

One Lyapunov function advocated for classic Kuramoto 
oscillators ( [2T] ) is the function W : Arc n (7r) — >• R defined 
for angles 6 in an open semi-circle and given by [52], [53] 



W(6) 



1 



Bi6 



(30) 



where B c G ]R nX ( n ( n - 1 )/ 2 ) i s an incidence matrix of the 
complete graph. As shown in [64, Theorem 4.4], the Lya- 
punov function ( [30] ) generalizes also to the coupled oscillator 
model (p}. Indeed, an even more general model is considered 
in [64], and a Lyapunov analysis yields the following result. 



Theorem 4.6: (Frequency synchronization I) Consider 
the coupled oscillator model (p} with a connected graph 
G(V,£,A) and uo G 1^. Assume that the algebraic con- 
nectivity is larger than a critical value, that is, 



\2(L) > Acriticai — ll^c 



(31) 



where B c G ]R nxn ( n - 1 )/ 2 i s the incidence matrix of the 
complete graph. Accordingly, define 7 max G ]7r/2,7r] and 
7min £ [0, 7r/2[ as unique solutions to (tt/2) • sinc(7 max ) = 
sin(7 m i n ) = A cr iticai/ A 2 (£). The following statements hold: 

1) phase cohesiveness: the set {0 G Arc n (7r) : 
|| B^Q || 2 < 7} C Ac (7) is positively invariant for ev- 
ery 7 G [7min, 7max], and each trajectory starting in the 
set {# G Arc n (7r) : ||i?j0|| < 7 ma x} asymptotically 
reaches the set {0 G Arc n (7r) : ||i?j0||2 < 7min}; and 

2) frequency synchronization: for every 0(0) G Arc n (7r) 
with 1 1 ^^6^(0) 1 1 2 < 7 m ax the frequencies 0i(t) synchro- 
nize exponentially to the average frequency cj sync = 
n^r=i^' anc *> gi yen phase cohesiveness in Ac (7) 
for some fixed 7 < tt/2, the exponential synchroniza- 
tion rate is no worse than Af e = —\2(L) cos(7). 

The proof of Theorem |4.6| follows a similar ultimate- 
boundedness strategy as the proof of Theorem |4.5| by using 
the Lyapunov function ( [30] ). It can be found in Appendix [B] 

For classic Kuramoto oscillators ( |2T] ), condition ( |3T] ) re- 
duces to K > || B^uj 1 1 2 . Clearly, the condition K > 1 1 B J 00 1 1 2 
is more conservative than the bound ( [26] ) which reads as if > 
^|| 00 — ^max — ^min- One reason for this conservatism 
is that the analysis leading to condition ( |3T] ) requires all 
phase distances 1 0i — 0j | to be bounded, whereas according 



to Lemma 4.2 only pairwise phase distances \0i — 0j 



{i,j} G £, need to be bounded for stable synchronization. 
The following result exploits these weaker assumptions and 
states a sharper (but only local) synchronization condition. 

Theorem 4.7: (Frequency synchronization II) Consider 
the coupled oscillator model (p} with a connected graph 
G(V, £ , A) and uj G 1^. There exists a locally exponentially 
stable equilibrium manifold [0] G Ag(tt/2) if 

A 2 (L) > \\B T u\\ 2 . (32) 

Moreover, if condition ( [32] ) holds, then [0] is phase cohesive 

in {0 eT n : \\B T 0\\ 2 < 7min } C A G ( 7min ), where 7min G 
[0,tt/2[ satisfies sin( 7min ) = ||£ t cj|| 2 /A 2 (L). 

The strategy to prove Theorem |4.7| is inspired by the 
ingenuous analysis in [52, Section IIV.B]. It relies on the 



insight gained from Lemma 4.2 that any synchronization 
manifold [0] G Ag*(7t/2) is locally stable, and it formulates 
the existence of such a synchronization manifold as a fixed 
point problem. Here, we follow the basic proof strategy in 
[52], but we provide a more accurate result together with a 
self-contained proof which is reported in Appendix [C] 

V. Conclusions and Open Research Directions 

In this paper we introduced the reader to the coupled os- 
cillator model (p}, we reviewed several applications, we dis- 
cussed different synchronization notions, and we presented 



different analysis approaches to phase synchronization, phase 
balancing, and frequency synchronization. 

Despite the vast literature, the countless applications, 
and the numerous theoretic results on the synchronization 
properties of model (p}, many interesting and important 
problems are still open. In the following, we summarize 
limitations of the existing analysis approaches and present 
a few worthwhile directions for future research. 

First, in many applications the coupling between the 
oscillators is not purely sinusoidal. For instance, phase delays 
in neuroscience [13], time delays in sensor networks [37], 
or transfer conductances in power networks [63] lead to a 
"shifted coupling" of the form sin(0i — 0j — (fij), where 
cfij G [— 7r/2, 7r/2]. In this case and also for other "skewed" 
or "symmetry-breaking" coupling functions, many of the 
presented analysis schemes either fail or lead to overly 
conservative results. Another interesting class of oscillator 
networks are systems of pulse-coupled oscillators featuring 
hybrid dynamics: impulsive coupling at discrete time instants 
and uncoupled continuous dynamics otherwise. This class of 
oscillator networks displays a very interesting phenomenol- 
ogy. For instance, the behavior of identical oscillators cou- 
pled in a complete graph strongly depends on the curvature 
of the uncoupled dynamics [142]. Most of the results known 
for continuously-coupled oscillators still need to be extended 
to pulse-coupled oscillators with hybrid dynamics. 

Second, in many applications [12], [24], [34], [63], [67] 
the coupled oscillator dynamics are not given by a simple 
first-order phase model of the form ([T]). Rather, the dynamics 
are of higher order, or sometimes there is no readily available 
phase variable to describe the limit cycle attracting the cou- 
pled dynamics. The analysis of oscillator networks with more 
general oscillator dynamics is largely unexplored. Whereas 
advances have been made for the simple case of phase 
synchronization of linear or passive oscillator networks, the 
case of frequency synchronization of non-identical oscillators 
with higher-order dynamics is not well-studied. 

Third, despite the vast scientific interest the quest for 
sharp, concise, and closed-form synchronization conditions 
for arbitrary complex graphs has been so far in vain [7], 
[8], [44]-[46], [74]. As suggested by Lemma |47T| Lemma 



4.2 Theorem 4.5 and the proof of Theorem |4.7| the proper 



metric for the synchronization problem is the incremental 
oo-norm y.B^lloo = max{ i5j } G £ \0i — 0j\. In the authors' 
opinion, a Banach space analysis of the coupled oscil- 
lator model ([T]) with the incremental oc-norm will most 
likely deliver the sharpest possible conditions. However, 
such an analysis is very challenging for arbitrary natural 
frequencies uj G 1^ and connected and weighted coupling 
graphs G(V, £, A). Recent work [114] by the authors puts 
forth a novel algebraic condition for synchronization with 
a rigorous analysis for specific classes of graphs and with 
(only) a statistical validation for generic weighted graphs. 

Fourth and finally, a few interesting and open theoretical 
challenges include the following. First, most of the presented 
analysis approaches and conditions do not extend to time- 
varying or directed coupling graphs G(V, £ , A), and alterna- 



tive methods need to be developed. Second, most known es- 
timates on the region of attraction of a synchronized solution 
are conservative. The semi-circle estimates given in Theorem 



4.3 and Theorem 4.5 rely on convexity of Arc n (7r) and are 
overly conservative. We refer to [63], [112] for a set of 
interesting results and conjectures on the region of attraction. 
Third, the presented analysis approaches are restricted to 
synchronized equilibria inside the set Ag(tt/2). Other inter- 
esting equilibrium configurations outside Ag(tt/2) include 
splay state equilibria or frequency-synchronized equilibria 
with phases spread over an entire semi-circle. 

We sincerely hope that this tutorial article stimulates 
further exciting research on synchronization in coupled oscil- 
lators, both on the theoretical side as well as in the countless 
applications. 

Appendix 

A. Modeling of the spring-interconnected particles 

Consider the spring network in Figure [T] consisting of a 
group of particles constrained to rotate around a circle of 
unit radius. For simplicity, we assume that the particles are 
allowed to move freely on the circle and exchange their order 
without collisions. Each particle is characterized by its phase 
angle 6i G S 1 and frequency 6i G R, and its inertial and 
damping coefficients are Mi > and Di > 0. 

The external forces and torques acting on each particle 
are (i) a viscous damping force DiOi opposing the direction 
of motion, (ii) a non-conservative force uji G R along the 
direction of motion depicting a preferred natural rotation 
frequency, and (iii) an elastic restoring torque between in- 
teracting particles i and j coupled by an ideal elastic spring 
with stiffness a^ > and zero rest length. 

To compute the elastic torque between the particles, we 
parametrize the position of each particle i by the unit vector 
Pi = [cos(0») , sin(0i)] T G S 1 C R 2 . The elastic Hookean 
energy stored in the springs is the function E : T n — » R 
given up to an additive constant by 



E(») 



E w ,«T""-»»> 

^2 ^ p dij (1 - cos(0i) cos(Oj) - sin(0i) sin((9 j )) 

E," 



an (l — cos(0i — 0j)) , 

where we employed the trigonometric identity cos (a — 0) = 
cos a cos /3+sin a sin j3 in the last equality. Hence, we obtain 
the restoring torque acting on particle i as 

d 

TAQ) = - — E(0) = - V an sin(^ - 6> 7 ) . 

Therefore, the network of spring-interconnected particles 
depicted in Figure [T] obeys the dynamics 



MSi + DiOi = uji - V dij sin(0i - 0j) , 

zG{l,...,n}. (33) 

The coupled oscillator model ([T]) is then obtained as the 
kinematic variant or the overdamped limit of the spring 



network ( [33] ) with zero inertia Mi = and unit damping 
Di = 1 for all oscillators i G {1, . . . , n}. 

B. Proof of Theorem \4.6\ 

Assume that 0(0) G Arc n (p) for p G [0, 7r[. Recall that the 
angular differences are well defined for in the open semi- 
circle Arc n (7r), and define the vector of phase differences 
S 4 BJO = (0 2 - l5 . . . ) G [-7T, +^] n(n - 1)/2 . By taking 
the derivative d/dtS(t) the phase differences satisfy 

5 = BJuj - BjBdmg({a ZJ } {Zjj} es) sin(B T 6) 
= B^uj - B T C B C diag({a ii } ijiG{ i > ... >n } ji<:7 -) sin(^), (34) 

where sin(x) = (sin(xi), . . . , sin(x n )) for a vector x G R n . 
Notice that for 0(0) G Arc n (7r) the ^-dynamics ( [34] ) are well- 
defined for an open interval of time. In the following, we 
will show that the set {5 G R n : ||5||2 < 7max} is positively 
invariant under condition ( [3T] ). As a consequence, the set 
{#GR n : ||5||oo < 7max < 7r} is positively invariant as 
well, and the ^-coordinates are well defined for alH > 0. 

The Lyapunov function ( [30] ) reads in ^-coordinates as 
W(S) = ^||^|| 2 , and its derivative along trajectories of ([34]) is 

x'T tz>T , s T B^B c dmg({a ij }i <j )sin(5) 



Hj fi<j J 

n S T diag({a i:) - } i<7 ) sm(6) , (35) 



W(S) = 5 r B- to 
= 5 T B^lu 

where the second equality follows from the identity 

5 T BjB c = T B c BjB c = T (nI n -l nxn )B c = nd T B c = nS. 
For || £2 1| < p, p e [0, 7r[, consider the following inequalities 

n5 T diag({o ij } i<j ) sin(5) 

= n (Bj8) T dmg({ aij sinc(^ - ^)} i<i )(B c T ^) 
> nsmc(p) (B c T ^) T diag({a ii } i<i )(B c T ^) 
> A 2 (L) sinc(p) \\B T c 9\\l = A 2 (L) sinc^^Hl , 

where the last inequality follows from [64, Lemma 4.7]. 
Hence, the derivative ( [35] ) simplifies further to 

W(S) < S t BJuj - A 2 (i)sinc(p)||(J|||. (36) 

In the following we regard BJuj as external disturbance 
affecting the otherwise stable (5-dynamics ( [34] ) and apply 
ultimate boundedness arguments [139]. Note that the right- 
hand side of ( [36] ) is strictly negative for 

ii cm . A \\Bjuj\\ 2 ^ ^critical 

11 112 > ^ C A 2 (L)sinc(p) A 2 (L)sinc(p) ' 

Pick e G]0, 1[. If p > \\S\\ 2 > n c /e, then the right-hand side 
of ( [36] ) is upper-bounded by 

W{5) < - (1 - e) • A 2 (L) smc(p)W{5) . 

In the following, choose p such that p > p > p c and let e = 
Pel p G ]0, 1[. By standard ultimate boundedness arguments 
[139, Theorem 4.18], for p(0)|| 2 < p, there is T > such 
that || 5(t) || 2 is exponentially decaying for t G [0, T] and 
\\S(t)\\ 2 < p for all t > T. For the choice p = 7 with 
7 G [0,7r/2[, the condition p > p c reduces to 



7sinc(p) > A crit i ca i/A 2 (L) . 



(37) 



Now, we perform a final analysis of the bound ( [37] ). The 
left-hand side of ( [37] ) is an increasing function of 7 and 
a decreasing function of p. Therefore, there exists some 
(p, 7) in the convex set A = {(p, 7) : p G [0,7r[, 7 G 
[0, 7r/2[, p > 7} satisfying equation ( [37] ) if and only if 
the inequality ( [37] ) is true atp = 7 = 7r/2, where the left- 
hand side of ( [37] ) achieves its supremum in A. The latter 
condition is equivalent to inequality ( [3T] ). Additionally, if 
these two equivalent statements are true, then there is an 
open set of points in A satisfying ( [37] ), which is bounded 
by the unique curve that satisfies inequality ( [37] ) with the 
equality sign, namely f(p, 7) = 0, where / : A —> R, 
/(P»7) = 7 sinc(p)- A critiC ai/ A 2 (i). Consequently, for every 
G {(p, 7 ) G A : /(p, 7 ) > 0}, it follows for 
\\5(0) || 2 < P that there is T > such that ||(S(t)|| 2 < 7 for all 
t >T. The supremum value for p is given by p max G ]7r/2, tt] 
solving the equation /(p ma x, tt/2) =0 and the infimum value 
of 7 by 7min e [0, tt/2[ solving the equation /(7 min , 7min) = 0. 

This proves statement 1) (where we replaced p max by 7 max ) 
and shows that there is T > such that WB^O^W^ < 
\\Bje(t)\\ 2 < 7min < tt/2 for all t > T. Thus, 0(t) G 
^G(lmm) for t >T, and frequency synchronization can be 
established analogously to the proof of Theorem |4.5| 



C. Proof of Theorem \4. 7\ 

According to Lemma |4~2] there exists a locally exponen- 
tially stable synchronization manifold [6] G Aq (7), 7 G 
[0,7r/2[, if and only if there is an equilibrium 6 G Ag(j). 
The equilibrium equations ( [T6] ) can be rewritten as 



L(B T 6)6, 



(38) 



j)}{i,i}e£:)^ T is 



where L(B T 0) = B diag({a^ sinc(#i 
the Laplacian matrix associated with the graph G(V,£,A) 
with nonnegative edge weights = sinc(#i — #j) for 
# G Ac (7). Since for any weighted Laplacian matrix L, we 
have that L = • L = I n — (l/n)l nxn (follows from 
the singular value decomposition [117]), a multiplication of 
equation ([38]) from the left by B T L(B T 6)^ yields 



B t L(B t 6) ] uj = B t Q. 



(39) 



Note that the left-hand side of equation ( [39] ) is a continuou^] 
function for # G Ac (7). Consider the formal substitution 
x = B T 0, the compact and convex set £00(7) = {% G 
£? T R n : II^Hoo < 7}, and the continuous map / : £00(7) —> 
R given by f(x) = B T L(x)^uj. Then equation ( [39] ) reads 
as the fixed-point equation /(x) = x, and we can invoke 
Brouwers's Fixed Point Theorem which states that every 
continuous map from a compact and convex set to itself has 
a fixed point, see for instance [143, Section 7, Corollary 8]. 

Since the analysis of the map / in the oo-norm is very 
hard in the general case, we resort to a 2-norm analysis and 
restrict ourselves to the set £2(7) = { x £ B T R n : ||x||2 < 
7} ^ £00(7)- By Brouwer's Fixed Point Theorem, there 

4 The continuity can be established when re-writing equations {38) and 
( [39) in the quotient space 1^, where L(B T 6) is nonsingular, and using 
the fact that the inverse of a matrix is a continuous function of its elements. 



exists a solution x G £2(7) to the equations x = f(x) if and 
only if || f{x) \\ 2 < 7 for all x G £2(7), or equivalently if 
and only if 



max \\B T L(xYuj\\ n < j . 
xes 2 { 1 ) 11 112 



(40) 



In the following we show that ( |32j > is a sufficient condition 
for inequality ( |4Q| . 

First, we establish some identities. For a Laplacian matrix 
L, we obtain = Fdiag(0, {l/Xi(L)} i=2 ... n )V T , where 
Ai(L) = and \i(L) > 0, i G {2, . . . , n}, are the eigenval- 
ues of L and V G M n x n is an associated orthonormal matrix 
of eigenvectors. It follows that Fdiag (0, 1, . . . , 1) V T = 
I n — (l/n)l nxn , and since uj _L l n , there exists a G R^l 
(not necessarily unique), such that uj = Ba. By means of 
these identities, the left-hand side of ( [40] ) can be simplified 
and upper-bounded for all x G £2(7): 



B t L(x) ] uj\\ 
B T V(x) diag ( 



B T L(x) ] Ba\ 
1 



V T (x)Ba 



< 



1 



X 2 (L(x)) 
(1/A 2 (i(x))) 



A 2 (L(x))''"' A n (L(x)) 
5 T V(x) diag (0, 1, ... , 1) V T (x)Ba 



B'uj 



(41) 



Thus, a sufficient condition for inequality ( |40| ) to be true can 
be derived as follows: 



max ||£> T L(£)^co>| 



2 < \\B T u\\ 2 max (l/A 2 (L(x)))) 



< 



|B T . 



max 



(l/\ 2 (L(x)))) 



= ||5 T ^|| 2 /(A 2 (L)-smc(7)) < 7 , 

where we used identity ( |4T] ), we enlarged the domain £2(7) 
to{xGM |£:| : Halloo < 7}, and we used the fact X 2 (L(x)) > 
X 2 (L) • sinc(7) for H^Hoo < 7. In summary, we conclude 
that there is a locally exponentially stable synchronization 
manifold [0] G {0 G T n : ||£ T <9|| 2 <7} C A G ( 7 )if 



A 2 (L)sin( 7 ) > 



a; 2 . 



(42) 



Since the left-hand side of ( [42] ) is a concave function of 7 G 
[0,7r/2[, there exists an open set of 7 G [0, 7r/2[ satisfying 
equation ( [42] ) if and only if equation ( |42| ) is true with the 
strict equality sign at 7* = tt/2, which corresponds to 
condition ( [32] ). Additionally, if these two equivalent state- 
ments are true, then there exists a unique 7 m i n G [0, 7r/2[ 
that satisfies equation ( [27] ) with the equality sign, namely 
sin(7 m i n ) = ||£? T co>||2/A2(I/). This concludes the proof. 
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